ts9

jholiday <- function(year, holiday.names = TRUE){
  if(missing(year)) stop("year was not given")
  stopifnot(is.numeric(year))
  stopifnot(length(year) == 1)
  stopifnot(year > 1948)
  ## We need strings of weekdays for a different locale.
  ## Starts with Sunday
  wdayStrings <- weekdays(seq(as.Date("2013/2/17"), by = "days", length = 7))

  .fixedDate <- function(x, name){
    h <- as.Date(paste(year, x, sep = "-"))
    if(holiday.names) names(h) <- name
    return(c(d, h))
  }
  .unfixedDate <- function(h, name){
    if(holiday.names) names(h) <- name
    return(c(d, h))
  }
  d <- as.Date(vector())                # There are no Date()
  n <- character()
  ## ====== January ==========
  # New Year's Day
  d <- .fixedDate("01-01", "New Year's Day")
  # Coming of Age Day
  if(year >= 2000){
    d <- .unfixedDate(.findDateByWeekday(year, 1, wdayStrings[2], 2),
                      "Coming of Age Day")
  }else{
    d <- .fixedDate("01-15", "Coming of Age Day")
  }

  ## ====== February ==========
  # Foundation Day
  if(year >= 1967){
    d <- .fixedDate("02-11", "Foundation Day")
  }else if(year == 1989){
    d <- .fixedDate("02-24", "State Funeral of the Showa Emperor")
  }

  ## ====== March ==========
  # Vernal Equinox Day
  d <- .unfixedDate(as.Date(paste(year, "03", .Shunbun(year), sep = "-")),
                    "Vernal Equinox Day")

  ## ====== April ==========
  # Showa Day
  if(year >= 2007){
    d <- .fixedDate("04-29", "Showa Day")
  }else if(year >= 1989){
    d <- .fixedDate("04-29", "Greenery Day")
  }else {
    d <- .fixedDate("04-29", "The Emperor's Birthday")
  }
  # Marriage of Crown Prince Akihito
  if(year == 1959){
    d <- .fixedDate("04-10", "Marriage of Crown Prince Akihito")
  }
  ## ====== May ==========
  # Constitution Memorial Day
  d <- .fixedDate("05-03", "Constitution Memorial Day")
  if(year >= 2007){
    d <- .fixedDate("05-04", "Greenery Day")
  }else if(year >= 1988){
    d <- .fixedDate("05-04", "Citizens' Holiday")
  }
  # Children's Day
  d <- .fixedDate("05-05", "Children's Day")
  
  ## ====== June ==========
  if(year == 1993){
    d <- .fixedDate("06-09", "Marriage of Crown Prince Naruhito")
  }

  ## ====== July ==========
  # Marine Day
  if(year >= 2003){
    d <- .unfixedDate(.findDateByWeekday(year, 7, wdayStrings[2], 3),
                      "Marine Day")
  } else if(year >= 1996){
    d <- .fixedDate("07-20", "Marine Day")
  }

  ## ====== August ==========
  # Mountain Day
  if(year >= 2016){
    d <- .fixedDate("08-11", "Mountain Day")
  }
  
  ## ====== September ==========
  # Autumnal Equinox Day
  aed <- as.Date(paste(year, "09", .Shubun(year), sep = "-"))
  d <- .unfixedDate(aed, "Autumnal Equinox Day")
  # Respect-for-the-Aged Day
  if(year >= 2003){
    rad <- .findDateByWeekday(year, 9, wdayStrings[2], 3)
    d <- .unfixedDate(rad, "Respect-for-the-Aged Day")
    if((aed - rad) == 2){
      d <- .unfixedDate(aed - 1, "Citizens' Holiday")
    }
  }else if(year >= 1966){
    d <- .fixedDate("09-15", "Respect-for-the-Aged Day")
  }

  ## ====== October ==========
  # Health and Sports Day
  if(year >= 2000){
    d <- .unfixedDate(.findDateByWeekday(year, 10, wdayStrings[2], 2),
                      "Health and Sports Day")
  }else if(year > 1966){
    d <- .fixedDate("10-10", "Health and Sports Day")
  }
  
  ## ====== November ==========
  # Culture Day
  d <- .fixedDate("11-03", "Culture Day")
  # Labour Thanksgiving Day
  d <- .fixedDate("11-23", "Labour Thanksgiving Day")
  # Official Enthronement Ceremony of Emperor Akihito
  if(year == 1990){
    d <- .fixedDate("11-12", "Official Enthronement Ceremony of Emperor Akihito")
  }

  ## ====== December ==========
  # The Emperor's Birthday
  if(year >= 1989){
    d <- .fixedDate("12-23", "The Emperor's Birthday")
  }

  ## ====== Others ==========
  # Transfer Holiday
  sun <- weekdays(d) == wdayStrings[1]
  if(any(sun)){
    hm <- d[sun] + 1
    if(year >= 2007){
      while(any(i <- hm %in% d)){
        hm[i] <- hm[i] + 1
      }
      if(holiday.names) names(hm) <- rep("Transfer Holiday", length(hm))
      d <- c(d, hm)
    }else if(year >= 1973){
      hm <- hm[! hm %in% d]
      if(holiday.names) names(hm) <- rep("Transfer Holiday", length(hm))
      d <- c(d, hm)
    }
  }
  return(sort(d))
}

is.jholiday <- function(dates){
  if(missing(dates)) stop("dates was not given")
  stopifnot(class(dates) == "Date")
  sapply(dates, function(x){
    y <- as.numeric(format(x, "%Y"))
    h <- jholiday(y, holiday.names = FALSE)
    x %in% h})
}

### ====== Fix me =======
### improve calculation cost
### =====================
.findDateByWeekday <- function(year, month, weekday, ordinal){
  stopifnot(is.numeric(year))
  stopifnot(is.numeric(month))
  stopifnot(is.character(weekday))
  stopifnot(is.numeric(ordinal))
  from <- as.Date(paste(year, month, 1, sep="/"))
  to <- as.Date(paste(year, month + 1, 1, sep="/")) - 1
  d <- seq(from, to, by="days")
  w <- weekdays(d) == weekday
  x <- d[w]
  return(x[ordinal])
}
  
## Calculation formulas in .Shunbun() and .Shubun() were referred to
## Shinkoyomi-benricho, Ephemeris Computation Workshop eds, Koseisha
## Koseikaku: Tokyo, 1991, ISBN: 9784769907008.

.Shunbun <- function(year){
  if(year <= 1947){
    dd = NULL;
  }
  else if(year <= 1979){
    dd = trunc(20.8357 + 0.242194 * (year - 1980) - trunc((year - 1983) / 4))
  }
  else if(year <= 2099){
    dd = trunc(20.8431+ 0.242194 * (year - 1980) - trunc((year - 1980) / 4))
  }
  else if(year <= 2150){
    dd = trunc(21.851 + 0.242194 * (year-1980) - trunc((year - 1980) / 4))
  }
  else {
    dd = NULL
  }
  return(dd)
}

.Shubun <- function(year){
  if(year <= 1947){
    dd = NULL;
  }
  else if(year <= 1979){
    dd = trunc(23.2588 + 0.242194 * (year - 1980) - trunc((year - 1983) / 4))
  }
  else if(year <= 2099){
    dd = trunc(23.2488+ 0.242194 * (year - 1980) - trunc((year - 1980) / 4))
  }
  else if(year <= 2150){
    dd = trunc(24.2488 + 0.242194 * (year-1980) -  trunc((year - 1980) / 4))
  }
  else {
    dd = NULL
  }
  return(dd)
}

#線形・ガウス型状態空間モデルにおける代表的な成分モデルの紹介と分析例

##個別のモデルの組み合わせ

##ローカルレベルモデル

###例: 人工的なローカルレベルモデル

#【ローカルレベルモデルに従う人工的なデータの作成】
# 前処理
set.seed(23)
library(dlm)
# ローカルレベルモデルの設定
W <- 1
V <- 2
m0 <- 10
C0 <- 9
mod <- dlmModPoly(order = 1, dW = W, dV = V, m0 = m0, C0 = C0)
# カルマン予測を活用して観測値を作成
t_max <- 200
sim_data <- dlmForecast(mod = mod, nAhead = t_max, sampleNew = 1)
y <- sim_data$newObs[[1]]
# 結果をts型に変換
y <- ts(as.vector(y))
# 結果のプロット
plot(y, ylab = "y")

##ローカルトレンドモデル

##周期モデル

###時間領域からのアプローチ

###周波数領域からのアプローチ

###例: 大気中の二酸化炭素濃度

####ローカルトレンドモデル+周期モデル(時間領域アプローチ)

#【ローカルトレンドモデル+周期モデル(時間領域アプローチ)】
# 前処理
library(dlm)
# データの読み込み
Ryori <- read.csv("https://raw.githubusercontent.com/hagijyun/tsbook/master/CO2.csv")
# データをts型に変換し、2014年12月までで打ち切る
y_all <- ts(data = Ryori$CO2, start = c(1987, 1), frequency = 12)
y <- window(y_all, end = c(2014, 12))
# モデルの設定：ローカルトレンドモデル+周期モデル(時間領域アプローチ)
build_dlm_CO2a <- function(par) {
  return(
    dlmModPoly(order = 2, dW = exp(par[1:2]), dV = exp(par[3])) +
    dlmModSeas(frequency = 12, dW = c(exp(par[4]), rep(0, times = 10)), dV = 0)
  )
}
# パラメータの最尤推定と結果の確認
fit_dlm_CO2a <- dlmMLE(y = y, parm = rep(0, 4), build = build_dlm_CO2a)
fit_dlm_CO2a
## $par
## [1]  -2.50345183 -21.91152657  -0.09834039  -5.18274450
## 
## $value
## [1] 330.4628
## 
## $counts
## function gradient 
##       39       39 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod  <- build_dlm_CO2a(fit_dlm_CO2a$par)
# カルマンフィルタリング
dlmFiltered_obj  <- dlmFilter(y = y, mod = mod)
dlmFiltered_obja <- dlmFiltered_obj              # 後で予測値を比較するために別名で保存
# フィルタリング分布の平均
   mu <- dropFirst(dlmFiltered_obj$m[, 1])
gamma <- dropFirst(dlmFiltered_obj$m[, 3])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値")
ts.plot(   mu, ylab = "レベル成分", ylim = c(350, 405))
ts.plot(gamma, ylab = "周期成分"  , ylim = c( -9,   6))
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] -330.4628

####ローカルレベルモデル+周期モデル(時間領域アプローチ)

#【ローカルレベルモデル+周期モデル(時間領域アプローチ)】
# モデルの設定：ローカルレベルモデル+周期モデル(時間領域アプローチ)
build_dlm_CO2b <- function(par) {
  return(
    dlmModPoly(order = 1, dW = exp(par[1]), dV = exp(par[2])) +
    dlmModSeas(frequency = 12, dW = c(exp(par[3]), rep(0, times = 10)), dV = 0)
  )
}
# 以降のコードは表示を省略
# パラメータの最尤推定と結果の確認
fit_dlm_CO2b <- dlmMLE(y = y, parm = rep(0, 3), build = build_dlm_CO2b)
fit_dlm_CO2b
## $par
## [1] -1.202149 -0.329150 -4.697639
## 
## $value
## [1] 340.1425
## 
## $counts
## function gradient 
##       16       16 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod  <- build_dlm_CO2b(fit_dlm_CO2b$par)
# カルマンフィルタリング
dlmFiltered_obj  <- dlmFilter(y = y, mod = mod)
dlmFiltered_objb <- dlmFiltered_obj              # 後で予測値を比較するために別名で保存
# フィルタリング分布の平均
   mu <- dropFirst(dlmFiltered_obj$m[, 1])
gamma <- dropFirst(dlmFiltered_obj$m[, 2])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値")
ts.plot(   mu, ylab = "レベル成分", ylim = c(350, 405))
ts.plot(gamma, ylab = "周期成分"  , ylim = c( -9,   6))
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] -340.1425

####ローカルトレンドモデル+周期モデル(周波数領域アプローチ)

#【ローカルトレンドモデル+周期モデル(周波数領域アプローチ)】
# モデルの設定：ローカルトレンドモデル+周期モデル(周波数領域アプローチ)
build_dlm_CO2c <- function(par) {
  return(
    dlmModPoly(order = 2, dW = exp(par[1:2]), dV = exp(par[3])) +
    dlmModTrig(s = 12, q = 2, dW = exp(par[4]), dV = 0)
  )
}
# 以降のコードは表示を省略
# パラメータの最尤推定と結果の確認
fit_dlm_CO2c <- dlmMLE(y = y, parm = rep(0, 4), build = build_dlm_CO2c)
fit_dlm_CO2c
## $par
## [1]  -2.75883291 -18.23154879   0.07782582  -7.96886928
## 
## $value
## [1] 284.6092
## 
## $counts
## function gradient 
##       29       29 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod  <- build_dlm_CO2c(fit_dlm_CO2c$par)
# カルマンフィルタリング
dlmFiltered_obj  <- dlmFilter(y = y, mod = mod)
dlmFiltered_objc <- dlmFiltered_obj              # 後で予測値を比較するために別名で保存
# フィルタリング分布の平均
   mu <- dropFirst(dlmFiltered_obj$m[, 1])
gamma <- dropFirst(dlmFiltered_obj$m[, 3] + dlmFiltered_obj$m[, 5])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値")
ts.plot(   mu, ylab = "レベル成分", ylim = c(350, 405))
ts.plot(gamma, ylab = "周期成分"  , ylim = c( -9,   6))
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] -284.6092

#【2015年からの予測】
# カルマン予測
dlmForecasted_object <- dlmForecast(mod = dlmFiltered_obj, nAhead = 12)
# 予測値の標準偏差・2.5%値・97.5%値を求める
f_sd <- sqrt(as.numeric(dlmForecasted_object$Q))
f_lower <- dlmForecasted_object$f + qnorm(0.025, sd = f_sd)
f_upper <- dlmForecasted_object$f + qnorm(0.975, sd = f_sd)
# 全観測値、予測値の平均値・2.5%値・97.5%値をts型として結合する
y_union <- ts.union(y_all, dlmForecasted_object$f, f_lower, f_upper)
# 以降のコードは表示を省略
# 結果のプロット
plot(y_union, plot.type = "single",
     xlim = c(2010, 2016),
     ylim = c( 385,  410), ylab = "", 
     lty = c("solid", "solid", "dashed", "dashed"),
     col = c("lightgray", "black", "black", "black"))
# 凡例
legend(legend = c("観測値", "平均 (予測分布)", "95%区間 (予測分布)"),
       lty = c("solid", "solid", "dashed"),
       col = c("lightgray", "black", "black"),
       x = "topleft", cex = 0.6)

#【2015年からの予測を3つのモデルで比較】
# モデルa, b, cの各々に対して、予測値の平均値・2.5%値・97.5%値を求める
f_all <- lapply(list(dlmFiltered_obja, dlmFiltered_objb, dlmFiltered_objc),
                function(x){
  # カルマン予測
  dlmForecasted_object <- dlmForecast(mod = x, nAhead = 12)
  # 予測値の標準偏差・2.5%値・97.5%値を求める
  f_sd <- sqrt(as.numeric(dlmForecasted_object$Q))
  f_lower <- dlmForecasted_object$f + qnorm(0.025, sd = f_sd)
  f_upper <- dlmForecasted_object$f + qnorm(0.975, sd = f_sd)
  
  # 結果をまとめる
  return(ts.union(
     mean = dlmForecasted_object$f,
    lower = f_lower,
    upper = f_upper
  ))
})
# 各モデルの予測結果をts型として統合する
names(f_all) <- c("a", "b", "c")
y_pred <- do.call("ts.union", f_all)
# 全観測値から2015年のデータを切り出す
y_true <- window(y_all, start = 2015)
# 以降のコードは表示を省略
# 結果のプロット
plot(y_pred, plot.type = "single", type = "b",
     xlab = "Time (2015年)", xaxt = "n", ylab = "",
     pch = c(rep("a", 3), rep("b", 3), rep("c", 3)),
     lty = rep(c("solid", "dashed", "dashed"), 3),
     col = rep(c("lightgray", "darkgray", "darkgray"), 3))
lines(y_true)
axis(side = 1, at = time(y_true), labels = 1:12)
# 凡例
legend(legend = c("観測値", "平均 (予測分布)", "95%区間 (予測分布)"),
       lty = c("solid", "solid", "dashed"),
       col = c("black", "lightgray", "darkgray"),
       x = "bottomleft", cex = 0.6)

#【2015年からのMAPEを3つのモデルで比較】
MAPE(true = y_true, pred = y_pred[, "a.mean"])
## [1] 0.001990843
MAPE(true = y_true, pred = y_pred[, "b.mean"])
## [1] 0.002313552
MAPE(true = y_true, pred = y_pred[, "c.mean"])
## [1] 0.001923007

##ARMAモデル

###例: 国産ビールの生産高

#【国産ビールの生産高】
# 前処理
library(dlm)
# データの読み込み
beer <- read.csv("https://raw.githubusercontent.com/hagijyun/tsbook/master/BEER.csv")
# データをts型に変換
y <- ts(beer$Shipping_Volume, frequency = 12, start = c(2003, 1))
# プロット
plot(y)

# データの対数変換
y <- log(y)
# 対数変換後のデータのプロット
plot(y, ylab = "log(y)")

####ローカルトレンドモデル+周期モデル(時間領域アプローチ)

#【国産ビールの生産高：ローカルトレンドモデル+周期モデル(時間領域アプローチ)】
# モデルの設定：ローカルトレンドモデル+周期モデル(時間領域アプローチ)
build_dlm_BEERa <- function(par){
  return(
    dlmModPoly(order = 2, dW = exp(par[1:2]), dV = exp(par[3])) +
    dlmModSeas(frequency = 12, dW = c(exp(par[4]), rep(0, times = 10)), dV = 0)
  )
}
# パラメータの最尤推定と結果の確認
fit_dlm_BEERa <- dlmMLE(y = y, parm = rep(0, 4), build = build_dlm_BEERa)
fit_dlm_BEERa
## $par
## [1] -10.519603 -18.095470  -5.645486  -9.843282
## 
## $value
## [1] -142.4162
## 
## $counts
## function gradient 
##       52       52 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod <- build_dlm_BEERa(fit_dlm_BEERa$par)
# カルマン平滑化
dlmSmoothed_obj <- dlmSmooth(y = y, mod = mod)
# 平滑化分布の平均
   mu <- dropFirst(dlmSmoothed_obj$s[, 1])
gamma <- dropFirst(dlmSmoothed_obj$s[, 3])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値(対数変換後)")
ts.plot(   mu, ylab = "レベル成分")
ts.plot(gamma, ylab = "周期成分")
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] 142.4162

####ローカルレベルモデル+周期モデル(時間領域アプローチ)

#【国産ビールの生産高：ローカルレベルモデル+周期モデル(時間領域アプローチ)】
# モデルの設定：ローカルレベルモデル+周期モデル(時間領域アプローチ)
build_dlm_BEERb <- function(par){
  return(
    dlmModPoly(order = 1, dW = exp(par[1]), dV = exp(par[2])) +
    dlmModSeas(frequency = 12, dW = c(exp(par[3]), rep(0, times = 10)), dV = 0)
  )
}
# 以降のコードは表示を省略
# パラメータの最尤推定と結果の確認
fit_dlm_BEERb <- dlmMLE(y = y, parm = rep(0, 3), build = build_dlm_BEERb)
## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular
fit_dlm_BEERb
## $par
## [1] -8.731695 -5.665727 -9.770276
## 
## $value
## [1] -151.1444
## 
## $counts
## function gradient 
##       21       21 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod <- build_dlm_BEERb(fit_dlm_BEERb$par)
# カルマン平滑化
dlmSmoothed_obj <- dlmSmooth(y = y, mod = mod)
# 平滑化分布の平均
   mu <- dropFirst(dlmSmoothed_obj$s[, 1])
gamma <- dropFirst(dlmSmoothed_obj$s[, 2])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値(対数変換後)")
ts.plot(   mu, ylab = "レベル成分")
ts.plot(gamma, ylab = "周期成分")
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] 151.1444

####ローカルレベルモデル+周期モデル(時間領域アプローチ)+ARMAモデル

#【国産ビールの生産高：AR(1)成分の考慮】
# モデルの設定：ローカルレベルモデル+周期モデル(時間領域アプローチ)+AR(1)モデル
build_dlm_BEERc <- function(par){
  return(
    dlmModPoly(order = 1, dW = exp(par[1]), dV = exp(par[2]))           +
    dlmModSeas(frequency = 12, dW = c(exp(par[3]), rep(0, 10)), dV = 0) +
    dlmModARMA(ar = ARtransPars(par[4]), sigma2 = exp(par[5]))
  )
}
# 以降のコードは表示を省略
# パラメータの最尤推定と結果の確認
fit_dlm_BEERc <- dlmMLE(y = y, parm = rep(0, 5), build = build_dlm_BEERc)
## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular

## Warning in dlmLL(y = y, mod = mod, debug = debug): a numerically singular 'V' has
## been slightly perturbed to make it nonsingular
fit_dlm_BEERc
## $par
## [1] -8.920552e+00 -6.058190e+00 -9.491816e+00  4.901214e-05 -6.937754e+00
## 
## $value
## [1] -152.6346
## 
## $counts
## function gradient 
##       65       65 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
ARtransPars(fit_dlm_BEERc$par[4])
## [1] 4.901214e-05
# パラメータの最尤推定結果をモデルに指定
mod <- build_dlm_BEERc(fit_dlm_BEERc$par)
# カルマン平滑化
dlmSmoothed_obj <- dlmSmooth(y = y, mod = mod)
# 平滑化分布の平均
   mu <- dropFirst(dlmSmoothed_obj$s[,  1])
gamma <- dropFirst(dlmSmoothed_obj$s[,  2])
 arma <- dropFirst(dlmSmoothed_obj$s[, 13])
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(4, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値(対数変換後)")
ts.plot(   mu, ylab = "レベル成分")
ts.plot(gamma, ylab = "周期成分")
ts.plot( arma, ylab = "AR(1)成分")
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)
# 対数尤度の確認
-dlmLL(y = y, mod = mod)
## [1] 152.6346

##回帰モデル

###例: 任天堂の株価

#【任天堂の株価】
# 前処理
library(dlm)
# データの読み込み
NINTENDO <- read.csv("https://raw.githubusercontent.com/hagijyun/tsbook/master/NINTENDO.csv")
NINTENDO$Date <- as.Date(NINTENDO$Date)
NIKKEI225 <- read.csv("https://raw.githubusercontent.com/hagijyun/tsbook/master/NIKKEI225.csv")
NIKKEI225$Date <- as.Date(NIKKEI225$Date)
# 観測値と説明変数を設定
y      <- NINTENDO$Close
x_dash <- NIKKEI225$Close
# 以降のコードは表示を省略
# プロット
plot(x = NINTENDO$Date , y = y     , xlab = ""    , ylab = "",
     ylim = c(10695, 28220), type = "l", col = "lightgray")
par(new=T)
plot(x = NIKKEI225$Date, y = x_dash, xlab = "Time", ylab = "",
     ylim = c(10695, 28220), type = "l", lty = "dashed"   )
# 凡例
legend(legend = c("SP of NINTENDO", "Nikkei"),
       lty = c("solid", "dashed"),
       col = c("lightgray", "black"),
       x = "topleft", cex = 0.6)

#【任天堂の株価のベータ値】
# モデルの設定：回帰モデル
build_dlm_REG <- function(par) {
  dlmModReg(X = x_dash, dW = exp(par[1:2]), dV = exp(par[3]))
}
# パラメータの最尤推定と結果の確認
fit_dlm_REG <- dlmMLE(y = y, parm = rep(0, 3), build = build_dlm_REG)
fit_dlm_REG
## $par
## [1]  4.358726e-07 -5.186503e+00 -5.466001e-07
## 
## $value
## [1] 1233.212
## 
## $counts
## function gradient 
##        8        8 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# パラメータの最尤推定結果をモデルに指定
mod  <- build_dlm_REG(fit_dlm_REG$par)
str(mod)
## List of 11
##  $ m0 : num [1:2] 0 0
##  $ C0 : num [1:2, 1:2] 1e+07 0e+00 0e+00 1e+07
##  $ FF : num [1, 1:2] 1 1
##  $ V  : num [1, 1] 1
##  $ GG : num [1:2, 1:2] 1 0 0 1
##  $ W  : num [1:2, 1:2] 1 0 0 0.00559
##  $ JFF: num [1, 1:2] 0 1
##  $ JV : NULL
##  $ JGG: NULL
##  $ JW : NULL
##  $ X  : num [1:160, 1] 14405 14562 14088 14202 14087 ...
##  - attr(*, "class")= chr "dlm"
# カルマン平滑化
dlmSmoothed_obj <- dlmSmooth(y = y, mod = mod)
# 以降のコードは表示を省略
# 平滑化分布の平均と標準偏差を求める
beta <- dropFirst(dlmSmoothed_obj$s[, 2])
beta_sdev <- sqrt(dropFirst(
               sapply(dlmSvd2var(dlmSmoothed_obj$U.S, dlmSmoothed_obj$D.S), function(x){
                 diag(x)[2]
               })
             ))
# 平滑化分布の50%区間のために、25%値と75%値を求める
beta_quant <- list(beta + qnorm(0.25, sd = beta_sdev),
                   beta + qnorm(0.75, sd = beta_sdev))
# 結果のプロット
 plot(x = NINTENDO$Date, y = beta, type="l", ylim = c(0.75, 2.0),
      xlab = "Time", ylab = "Beta est")
lines(x = NINTENDO$Date, y = beta_quant[[1]], lty = "dashed")
lines(x = NINTENDO$Date, y = beta_quant[[2]], lty = "dashed")
# 凡例
legend(legend = c("平均 (平滑化分布)", "50%区間 (平滑化分布)"),
       lty = c("solid", "dashed"),
       col = c("black", "black"),
       x = "topleft", cex = 0.6)
# 参考イベント
mtext("×", at = as.Date("2015/3/17"), side = 1, adj = 0.5, line = -0.5)
mtext("×", at = as.Date("2016/7/6" ), side = 1, adj = 0.5, line = -0.5)
mtext("×", at = as.Date("2016/7/22"), side = 1, adj = 0.5, line = -0.5)
mtext("2015/3/17", at = as.Date("2015/3/17"), side = 1, adj = 0, cex = 0.6)
mtext("2016/7/6" , at = as.Date("2016/7/6" ), side = 1, adj = 1, cex = 0.6)
mtext("2016/7/22", at = as.Date("2016/7/22"), side = 1, adj = 0, cex = 0.6)

###例: ナイル川の流量(1899年の急減を考慮)

#【ナイル川の流量データにローカルレベルモデル+回帰モデル（干渉変数）を適用】
# 前処理
set.seed(123)
library(dlm)
# ナイル川の流量データ
y <- Nile
t_max <- length(y)
# 説明変数（干渉変数）を設定
x_dash <- rep(0, t_max)                  # 初期値として全て0（ダムなし）
x_dash[which(1899 <= time(y))] <- 1      # 1899年以後は全て1（ダムあり）
# ローカルレベルモデル+回帰モデル（干渉変数）を構築する関数
build_dlm_DAM <- function(par) {
  return(
    dlmModPoly(order = 1, dV = exp(par[1]), dW = exp(par[2])) +
    dlmModReg(X = x_dash, addInt = FALSE, dW = exp(par[3]), dV = 0)
  )
}
# パラメータの最尤推定
fit_dlm_DAM <- dlmMLE(y = y, parm = rep(0, 3), build = build_dlm_DAM)
modtv <- build_dlm_DAM(fit_dlm_DAM$par)
# カルマン平滑化
dlmSmoothed_obj <- dlmSmooth(y = y, mod = modtv)
# 平滑化分布の平均と分散
stv <- dropFirst(dlmSmoothed_obj$s)
stv_var <- dlmSvd2var(dlmSmoothed_obj$U.S, dlmSmoothed_obj$D.S)
stv_var <- stv_var[-1]
# 推定量の平均
s <- stv[, 1] + x_dash * stv[, 2]                       # x_dashも考慮
# レベル推定量の50%区間（25%値と75%値を求める）
coeff <- cbind(1, x_dash)
s_sdev <- sqrt(sapply(seq_along(stv_var), function(ct){ # 共分散も考慮
            coeff[ct, ] %*% stv_var[[ct]] %*% t(coeff[ct, , drop = FALSE])
          }))           
s_quant <- list(s + qnorm(0.25, sd = s_sdev), s + qnorm(0.75, sd = s_sdev))
# 以降のコードは表示を省略
# プロット
ts.plot(cbind(y, s, do.call("cbind", s_quant)),
        lty=c("solid", "solid", "dashed", "dashed"),
        col=c("lightgray", "black", "black", "black"))
# 凡例
legend(legend = c("観測値", "平均", "50%区間"),
       lty = c("solid", "solid", "dashed"),
       col = c("lightgray", "black", "black"),
       x = "topright", cex = 0.6)

###例: 家計における食費(曜日によって異なる影響を考慮)

#【家計における消費支出（食料）】
# 前処理
library(dlm)
# データの読み込み
food <- read.csv("https://raw.githubusercontent.com/hagijyun/tsbook/master/FOOD.csv")
# データをts型に変換
y <- ts(food$Expenditure, frequency = 12, start = c(2000, 1))
# プロット
plot(y)

# データの対数変換
y <- log(y)
# 対数変換後のデータのプロット
plot(y, ylab = "log(y)")

#【説明変数（平日効果）の設定】
# 日本の平休日を返すユーザ定義関数
jholidays <- function(days){
  # is.jholiday()を利用する
  library(Nippon)
  # 曜日（Day Of the Week）を求める
  DOW <- weekdays(days)
  # 土日祝（代休含む）は休日とみなす
  holidays <- (DOW %in% c("Saturday", "Sunday")) | is.jholiday(days)
  # 曜日を"休日" or "平日"で上書き
  DOW[ holidays] <- "Holidays"
  DOW[!holidays] <- "Weekdays"
  return(DOW)
}
# 検討対象期間の日付の数列
days <- seq(from = as.Date("2000/1/1"), to = as.Date("2009/12/31"), by = "day")
# 月毎に平休日の日数を集計
monthly <- table(substr(days, start = 1, stop = 7), jholidays(days))  
## Loading required package: stringr
## 
## Attaching package: 'Nippon'
## The following objects are masked _by_ '.GlobalEnv':
## 
##     is.jholiday, jholiday
# 説明変数（ある月の平日から休日の総数を引いた値）
x_dash_weekday <- monthly[, "Weekdays"] - monthly[, "Holidays"]

#【説明変数（うるう年効果）の設定】
# データ長
t_max <- length(y)
# 検討期間中のうるう年の2月
LEAPYEAR_FEB <- (c(2000, 2004, 2008) - 2000)*12 + 2
# 説明変数（うるう年の2月のみ1）
x_dash_leapyear <- rep(0, t_max)          # 初期値は全て0
x_dash_leapyear[LEAPYEAR_FEB] <- 1        # うるう年の2月は1

#【ローカルレベルモデル+周期モデル(時間領域アプローチ)+回帰モデルによる食費の分析】
# 説明変数（平日効果、うるう年効果）を統合
x_dash <- cbind(x_dash_weekday, x_dash_leapyear)
# ローカルレベルモデル+周期モデル（時間領域）+回帰モデルを構築する関数
build_dlm_FOODb <- function(par) {
  return(
    dlmModPoly(order = 1, dW = exp(par[1]), dV = exp(par[2]))         +
    dlmModSeas(frequency = 12, dW = c(0, rep(0, times = 10)), dV = 0) +
    dlmModReg(X = x_dash, addInt = FALSE, dV = 0)
  )
}
# パラメータの最尤推定
fit_dlm_FOODb <- dlmMLE(y = y, parm = rep(0, 2), build = build_dlm_FOODb)
# パラメータの最尤推定結果をモデルに指定
mod  <- build_dlm_FOODb(fit_dlm_FOODb$par)
-dlmLL(y = y, mod = mod)
## [1] 319.8928
# カルマンフィルタリング
dlmSmoothed_obj  <- dlmSmooth(y = y, mod = mod)
# 平滑化分布の平均
    mu <- dropFirst(dlmSmoothed_obj$s[, 1])
 gamma <- dropFirst(dlmSmoothed_obj$s[, 3])
beta_w <- dropFirst(dlmSmoothed_obj$s[, 13])[t_max]  # 時不変
beta_l <- dropFirst(dlmSmoothed_obj$s[, 14])[t_max]  # 時不変
# 結果の確認
cat(beta_w, beta_l, "\n")
## -0.003098287 0.03920194
# 回帰成分の平均
reg <- x_dash %*% c(beta_w, beta_l)
tsp(reg) <- tsp(y)
# 結果のプロット
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(4, 1)); par(oma = c(2, 0, 0, 0)); par(mar = c(2, 4, 1, 1))
ts.plot(    y, ylab = "観測値(対数変換後)")
ts.plot(   mu, ylab = "レベル成分")
ts.plot(gamma, ylab = "周期成分")
ts.plot(  reg, ylab = "回帰成分")
mtext(text = "Time", side = 1, line = 1, outer = TRUE)

par(oldpar)